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Abstract 

We apply the holonomic gradient method introduced by Nakayama et al. [3] 
to the evaluation of the exact distribution function of the largest root of a Wishart 
matrix, which involves a hypergeometric function \F\ of a matrix argument. Nu- 
merical evaluation of the hypergeometric function has been one of the longstanding 
problems in multivariate distribution theory. The holonomic gradient method offers 
a totally new approach, which is complementary to the infinite series expansion 
around the origin in terms of zonal polynomials. It allows us to move away from 
the origin by the use of partial differential equations satisfied by the hypergeomet- 
ric function. From numerical viewpoint we show that the method works well up 
to dimension 10. From theoretical viewpoint the method offers many challenging 
problems both to statistics and D-module theory. 

Keywords and phrases: -D-modules, Grobner basis, hypergeometric function of a matrix 
argument, zonal polynomial 

1 Introduction 



For multivariate distribution theory in statistics, the theory of zonal polynomials and 
hypergeometric functions of matrix arguments, introduced by A.T. James and other au- 
thors, was a very important development in the 1950's. They allowed explicit expressions 
of density functions and cumulative distribution functions of basic test statistics under 
non-null cases. Zonal polynomials are based on the representation theory of real general 
linear group and they possess many interesting combinatorial properties. Properties and 
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applications of zonal polynomials and hypergeometric functions of matrix arguments are 



surveyed in Gross and Richards [4J and Richards 21|. Zonal polynomials are special cases 



of Jack polynomials, whose properties have been intensively studied by many mathemati- 



cians. See for example Chapters VI and VII of Macdonald uM and Stanley |25j. Jack 



polynomials are further generalized to Macdonald polynomials (see, e.g., Kuznetsov and 
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Despite the above nice mathematical properties of zonal polynomials and hyperge- 
ometric functions of matrix arguments, from practical viewpoint they were not really 
useful for computations. Coefficients of zonal polynomials can be computed only through 
nontrivial combinatorial recursions. Although very ingenious recursion algorithms have 
been recently developed (Koev and Edelman 10] ), computing zonal polynomials of large 
degrees remains to be a difficult problem because of inherent combinatorial complexities. 
Also, the convergence of infinite series expansion of hypergeometric functions of a matrix 
argument in terms of zonal polynomials was found to be slow (Muirhead [l7j , Hashiguchi 
and Niki jsj]). Since the expansion of the hypergeometric function in terms of zonal poly- 
nomials is the expansion at the origin, the convergence for large values of the argument 
is necessarily slow. 

The holonomic gradient method allows us to move away from the origin by the use 
of partial differential equations. Thus our approach provides a promising new method 
for attacking a longstanding problem in multivariate statistics. Our holonomic gradient 
method is, in spirit, on the track of the holonomic systems approach to combinatorial 



identities by Zeilberger |32j. Note that the series expansion and our holonomic gradi- 
ent method are in fact complementary methods, because our method needs the series 
expansion for obtaining initial values for the partial differential equations. 

The main purpose of this paper is to verify the performance of holonomic gradient 
method for -J?\. We found that a straightforward implementation of the holonomic gra- 
dient method works well for dimensions up to 10. 

Butler and Wood ^ showed that the Laplace method gives a very good approximation 
to \F\ even for a high dimension, e.g., m = 32. However the Laplace method needs a 
peaked density function, which corresponds to a large degrees of freedom. Our method is 
an exact method, where the errors only come from discretization in numerically solving 
differential equations and the accuracies in the initial values. Hence our method works 
even for small degrees of freedom. 

The organization of this paper is as follows. In Section [2] we summarize preliminary 
facts on the exact distribution of the largest root of a Wishart matrix. In particular we 



state the partial differential equation for ±F 1 by Muirhead 16j. In Section[31 for expository 



purposes, we fully describe our holonomic gradient method for dimension two. In Section 
H] we derive properties of Pfaffian system for general dimensions. The Pfaffian system is 
a system of partial differential equations and is called an integrable connection in some 
literatures. Results of symbolic computations are presented in Section [5] and results of 
numerical experiments are presented in Section We end the paper with discussion of 
open problems in Section 
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2 Preliminaries 

Let k — (ki, . . . , ki) h k be a partition of a non- negative integer fc and define the Pochham- 
mer symbol (a) K by 

(a) K = nf a -^-j > (^k^Hia + j-l) ((a) = l). 

Let C K (y) denote the ( "C-normalization" of) zonal polynomial indexed by k of an m x m 
symmetric matrix Y . it is a homogeneous symmetric polynomial of degree k in the char- 
acteristic roots j/i, ... , 2/ m of y, satisfying V) hfe C K (y) = (trF) fc . For zonal polynomials 



in statistics see, e.g., James |8|], Muirhead [18j . Takemura |30| and Mathai et al. [15J. A 
hypergeometric function of a matrix argument is defined (Constantine jij) as 

Win n-r r . v ^ _ \^ \^ {^) n ■ ■ ■ {a p ) R C K {Y) ^ 

pT q\ a ly • • • j a p> C l> • • • i C qi r ) / / 



(ci) K . . . (c,)« fc! 



In this paper we study holonomic gradient method for iFi(a; c; Y). Let L m denote the 
m x m identity matrix and let |X| denote the determinant of X. For 3fta > (m + l)/2, 
9R(6 — a) > (m + l)/2, i-Fi(a; c; F) has the following integral representation 

xFM c; y) = - ["^ , / exp(trXy)|Xr-( m+1 )/ 2 |/ m - X\ c ~ a ~ W*dX, 

J- m{CL)L m {C — a) Jo<X<Im 

(2) 

where < X < I m means that X and I m — X are positive definite, dX = Yli<j m t ne 
Lebesgue measure of the upper triangular entries of X, and 



r m W=^I]r(a-Y) 

8=1 ^ ' 



The hypergeometric function xFi satisfies the the following Kummer relation (see (2.8) of 
Herz 0, (51) of James 0): 

exp(-try) 1 F 1 (a;c;y) = x F x (c- a,c; -Y). (3) 

Note that (j2J) implies that X F X is an entire function in Y. 

The cumulative distribution function of the largest root t x of the m x m Wishart 
matrix W with n degrees of freedom and the covariance matrix E is written as follows 



Pr[A<x] = Cexp(-| t rE-')>-» 1 F, (^±1; 2 + ^ + 1, £ E -') 



(4) 



where 

m \ 2 / 



This follows from the results in Section 9 of Constantine m and the Kummer relation fl3 



See also Sugiyama [27 . 

The following partial differential equations for iFi(a;b;Y) were derived by Muirhead 
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Theorem 1 (Theorem 5.1 of Muirhead [16j, Theorem 7.5.6 of Muirhead (18J). The hy- 

pergeometric function F = i-F\(a; c; Y) of a matrix argument Y = diag(yi, . . . , y rn ) is the 
unique solution of the following set of m partial differential equations 



m — 1 



1 m 



ID 



2 i=~* Vi ~ y > 



Vj 



2 i * - ^ 



5j — a 



F = 0, (5) 
(i = l, 



subject to the conditions that F is symmetric in yi, ■ ■ ■ ,y m an d F analytic at Y = 0, 
F(0) = 1. 

The partial differential equation (jSJ) has singularities along y,j = and yj = y^ j ^ i. 
However since F is an entire function, F is determined by the partial differential equations 



on the open region X = {y G ( 
the non-diagonal region. Using 



YYiLi Vi UijijiVi - Vi) 0}- In tnis paper we call X 



Vi 



1 + 



yi - Vj yi - yj 

we can rewrite ([5]) as giF — 0, i — 1, . . . , m, where 



i m 

9i = Vidi + (c - yi)^ + - J] ^ - dj) - a 



2 ^ - 



(6) 



is a differential operator annihilating F. In our holonomic gradient method we make a 
direct use of the partial differential equations for numerical evaluation of liq. 



3 Holonomic gradient method for dimension two 

In this section we illustrate the holonomic gradient method for the case of m = 2. Al- 
though our purpose is to implement an algorithm of our method for a larger dimension, 
for clarity it is best to do "by hand" calculation for the case of m = 2. As in the previous 
section we simply write F( Y) = iF(a; c; Y). 



In Nakayama et al. 19] the holonomic gradient method was used to obtain the max- 
imum likelihood estimate. The reciprocal of the likelihood function was minimized and 
the method was called the holonomic gradient descent. For the application of this paper 
we simply use the holonomic gradient method for evaluating F. Hence we omit the term 
"descent" . Also, for minimization, at each step of the iteration, a direction for increments 
was chosen to decrease the value of the function. In our application, starting from the 
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origin Y = 0, we can choose arbitrary path to the target value Y where we want to 
evaluate F(Y). 

Another minor difference of the expository explanation in this section from Nakayama 



et al. |19| and Sei et al. |24J is that we use the simple forward Euler method (e.g., Section 
3.1 of Ascher and Petzold [if) for updating partial derivatives of F. In Nakayama et al. 
19j , once an updating direction is chosen at each step of the iteration, the 4-th order 
Runge-Kutta method was used. The simple Euler method is used only for the purpose of 
exposition. It is easier to explain the basic idea of the holonomic gradient method with the 
simple Euler method. In our actual implementation in Section Owe use the Runge-Kutta 
method for numerically solving the differential equation. 

We will reduce our problem to a traditional problem of numerical analysis of an or- 
dinary differential equation (ODE). For the reduction we utilize the notion of holonomic 
differential equations and the gradients of their solutions. It is why we call our method 
holonomic gradient method. 

In the following we discuss the case of y\ ^ y 2 and y\ = y 2 separately. 

3.1 Holonomic gradient method for non-diagonal region 

In this subsection we assume y± ^ y 2 . Two partial differential equations in ([6]) are written 

as 

y x d\ + (c - y 1 )d 1 + (ft - ft) - al F = 0, (7) 

2 l/i - 2/2 J 



y 2 d\ + (c - y 2 )d 2 + — (ft - ft) - a 
2 2/2-2/1 



F = 0. 



Suppose that we want to evaluate a higher derivative d^d^F = d^d^F of F. Let 
n 2 > 2. Then by © 

d^d^F = d^d^~ 2 { - -ft + ft - \— ^ -Ad 2 - ft) + -)F (9) 

2/2 2y 2 (y 2 -y 1 ) y 2 ' 

Noting 

p. 1 1 rs 2/1 2/l(22/2-2/l) 

o 2 — — — O2 



'2/2 2/1' 2/2(2/2-2/1) 2/2(2/2-2/1; 
for n 2 > 2, the right-hand side of (j9j) is further written as 



^n 2 -sf^_ d2 _ c-Kq + 1 2/i(22/ 2 -2/i) 
-2/2 2/2 2y^(y 2 -y i y 



- \ - f yi M - ftft) - - 2 + -b 2 )f (io) 
22/2(2/2-2/1) 2/2 2/2 / 

Although the result is somewhat complicated, the important fact is that the total degree 
of differentiation n\ + n 2 on the left-hand side of (Q is decreased by one to n\ + n 2 — 1 
in fnU]) . As long as the degree of ft or ft is more than one, then we can recursively apply 
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(J7|) or (jSJ) to decrease the total degree of differentiation. It follows that for each ni,n 2 , 

(11) 



there exist rational functions h^' n2 \ h^' n2 \ h^' n2 \ h^ 1 ' 71 ^ in (2/1,2/2) such that 



dp&fF = h { ^ n2) F + hf^dxF + h^' n2) d 2 F + h^ 1 ' n2) d 1 d 2 F. 
In this notation (J7J) is written as 



dlF 



a_p _ / C-yi + 1 1/2 
2/i 2/1 22/1(2/! -y 2 ) 



)ftF + 



2/2 



2 2/i(yi - 2/2) 



ftF 



,(2,0) 



(2,0); 



(M 2 i 0) - o). 



(12) 



For a general dimension, (1111) corresponds to the reduction by a Grobner basis as discussed 
in Section HI 

For us the important case is n\ — 1, n 2 — 2. Since 



ft 



2/i 



2/2(2/2 - 2/1) 



4( 



— ! 1 )= 1 

2/2-2/1 2/2 (2/2-2/1) 2 



we have 
ftft 2 F = ft ( - 

= ( 



2/2 



ft 



2/2 

c- 2/2 
2/2 



22/2(2/2-2/1) 2/2 



Sift - - 



2(2/2-2/i) 2 



(ft - ft 



2/i 



,-(ftft-ft 2 ) + -ft)F 
22/2(2/2-2/1) 2/2 



There is a term 2/1 <9 2 on the right-hand side, into which we further substitute (171). Then 
( ITT]) for ftSfF is written as 



ftft 2 F 



2/2 



2/2 



Sift 



2 (2/2 - 2/i) 2 



(ft - ft) - I 



„ , rftft H Ci 

22/2(2/2-2/1) 2/2 



— 1 -((c - 2/i)ft + i^ 2 — (ft - ft) - a) )F 

22/2(2/2-2/1) 22/1-2/2 '/ 



22/2(2/2 - 2/i) 
3 1 



1 



c-2/1 



rftF 



4 (2/2 - 2/i) 2 2/2 22/2(2/2 - 2/i) 
c-2/2 , 1 2/i 



ftF 



4 (2/2 - 2/i) 2 "^ ^ 2/2 22/2(2/2-2/1) 
(1 ' 2) F + tfn^ftF + /^ 2) ftF + tfi 2) ftftF 



ftftF 



'on 



(13) 



Since F is a symmetric function in 2/1 and 2/2, dfd 2 F is obtained by permuting 2/1 and y 2 . 
Let 

/ ^ \ 

ftF 

ftF 
VftftF/ 
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denote the vector consisting of F and its square-free mixed derivatives. Differentiate the 
components of F by y\ and denote d\F = (diF, d 2 F, d\d 2 F, d\ d 2 Fy. Similarly define d 2 F. 
Then by (fT2l and (fT3"|). d { F, i = 1, 2, are written as d t F = Pi(Y)F, where Pi and P 2 are 
the following 4x4 matrices with rational function entries 



Pion 



( 


1 





o \ 


"oo 


"lO 


1.(2,0) 
""01 
















U (2 ' l} 

\"oo 


'no 


"oi 





P2(Y) 



( 





1 


^ 











1 


i(0,2) 
"00 


l(0,2) 

'no 


l(0,2) 
"01 





U (1,2) 

\"oo 


/i (1 ' 2) 
'no 


"oi 





The matrices P\,P 2 are called coefficient matrices of a Pfaffian system (an integrable 
connection) for F (Nakayama et al. [H|). Note that P 2 is obtained from Pi by permutation 
of yi and y 2 . If we know the values of the components of F at Y — (j/i, 2/2) 5 2/i 7^ 2/2j then 
values at a nearby point F + AY = (y 1 + Ay^y 2 + Ay 2 ) can be approximated by the 
simple Euler method (i.e. linear approximation) as 



F(Y + AY) = F{Y) + Ay 1 d 1 F{Y) + Ay 2 d 2 F(Y) 

= FiY) + A yi P x {Y)F{Y) + Ay 2 P 2 (Y)F(Y). 



(14) 



Now suppose that we want to evaluate F(yi,y 2 ) at a particular point (2/1,2/2) with 



2/i 7^ 2/2- If we know F(F ) at some point F = (yj 



(0) „(0), 
2/2 



(0) 



then we can choose an appropriate sequence of points = (yf\y2^), I 



close to the origin, 
0, . . . , L, such 

that (2/1,2/2) = {y'i J \y 2 J> )- Along the sequence we can use (fT4"|) to update F(F^) and 
finally the first element of P(F^ L ^) gives F(yi,y 2 ). 

Therefore it remains to consider how to obtain the initial values. Close to the origin 
we can use the definition (TjQ) of iF\. If F is very close to zero, then we only need zonal 
polynomials of low orders, whose explicit forms are known. Zonal polynomials up to the 
third order are as follows; Cm(Y) = A^(i)(F), 



'C(2)(F)\ 
^(1,1) (Y)J ~ 



(1 


-\ 




3 


( M {2) (Y)\ 


v° 


4 

3) 


\M {1>1) (Y)) 



C( 3 )(F) 
C(2,i)(Y) ] = 

vc ( i,i,i)(n 



/1 2 2 -\ 

5 5 

12 18 

T T 

\0 2 / 



A^(3)(F) 

^(14,1)00, 



(15) 

where Ai K (Y) is the monomial symmetric polynomial associated with a partition k. Since 
F(yi,y 2 ) can be expanded as 

(a)(2) 



F(2/i,2/ 2 ) = l 



( fl )(i) r m , 1 
— C(i)(F) + — 



C(2)(n + MM C(1 , 1)(n ) + ... 

( c J(i,i) / 



C)(2) ' (c)(l,l) 

l c )(l) 4 C )(2) \^l C J(2) ^l c J(l,l)/ 

(16) 
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for an example, <9i<9 2 F(0, 0) is obtained as 

d 1 d 2 F(0, 0) 
In a similar manner, we have 



(a) 2 2a(a - |) 



3(c) 2 3c(c-|)' 



9iF(0,0) = <9 2 F(0,0) = -, 9?F(0,0) = <9 2 F(0,0) 



[a) 2 



dld 2 F{0, 0) = <9 2 2 <9iF(0, 0) 



(«)3 



4(a) 2 (a- |) 



5(c) 3 5(c) 2 (c-| 



(17) 



These formulae can be obtained by a symbolic mathematics software, such as the routines 
for Jack polynomials in sage mathematics software system (Stein et al. (26[ 

In order to obtain the initial value F(Y ) at Y 
use the approximation 



(2/i°^ 5 2/2°^ ) dose to the origin, we can 



F(y?\yi 0) ] 



F(0, 0) + y[ o) d 1 F(0 : 0) + y 2 ">d 2 F(0, 0). 



(0); 



'11 



We code the above procedure using deSolve package in the data analysis system R. We 
show a simple source program in Appendix [B] In addition, since the zonal polynomials 
are easy to evaluate for m = 2, we also evaluate the series expansion of i F 1 up to k = 150. 
As an example, we compute percentage points by two methods for the case of n = 3, 
X = diag(l/2, 1/4). The following percentage points for £1 agree in two methods to 6 
digits. 



50% 



90% 



95% 



99% 



1.63785 3.54999 4.31600 6.05836 



3.2 Holonomic gradient method for the diagonal line 

In the previous subsection we assumed y± 7^ y 2 to avoid singularity of the differential 
equations. However ±F± itself does not have singularities. Hence we should be able to 
derive some differential equation even for y = yi — y 2 . 

In Q and (jSJ) we can perform the limiting operation y\ — > y 2 = y using the l'Hopital 
rule. Since F is a symmetric function, at (y, y) we have 

9iF(y,y) = d 2 F(y,y). 

Also dfF(y,y) = d 2 F(y,y). Hence by the l'Hopital rule, in (j7|) we have 



lim 

2/1-^3/2=3/ 



d x F - d 2 F 



yi - 2/2 

Then (J7|) for at (y, y) is written as 







y df+(c-y)d 1 + ^(dt-d 1 d 2 ) 



d\F - d x d 2 F = dfF - d^F. 



-yd 2 1 + {c-y)d 1 - y -d 1 d 2 



F. (19) 
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Based on this we derive an ODE for f(y) = F(y, y). Firstly, 

f\y) = 2d x F or d x F = f(y)/2. 



Secondly, 
From (|T9|) 

and 



or 



f"(y) = 2dfF + 2d x d 2 F 



\ydlF = ^yd 1 d 2 F - (c - y)d 1 F + aF, 



\yf\y) = \yd x d 2 F - (c - y)d x F + aF + ^yd x d 2 F 
= 2yd 1 d 2 F - (c - y)d 1 F + aF 
= 2yd 1 d 2 F- C —yf(y) + af(y) 



WFfay) = lf\y) + C -^-f(y) - Y y f{y) - (20) 



Thirdly, 



f"\y)= 2dlF + Zdld 2 F (21) 
In order to get another relation for d\F and dfd 2 F, we differentiate ([7]) by y 2 . Then by 

g 2 ^L_ = i) = - — — — — (22) 

2/1 -yi 2/1 - 2/2 (2/1 - 2/2r 

we obtain the following differential operator annihilating F: 

yid\d 2 + (c - yO^^ + 1 yi (fr - «9 2 ) + i^ 2 — (d^ - «9 2 ) - ad 2 . (23) 

2(2/i-2/ 2 ) 2 2y x -y 2 

Noting 2/2/(2/1 — 2/2) = 2/1/(2/1 — 2/2) — 1 this can be further written as 

v m + (c - m)BA + ^ (ft-ft) + ^-w)(fta-aS) _ 1 (aA _ ^ _ a92 

2 (2/1 - 2/2r 2 

= 2/i"i"2 + (c - 1 - 2/i)<9i<9 2 + — 7 V -{did 2 + <9 2 ) - ad 2 . 

2 (2/1 - 2/2) 2 2 

We now apply the l'Hopital rule to 

(d 1 -d 2 ) + (y 1 -y 2 )(d 1 d 2 -dl) 
(2/1 - 2/2) 2 
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We again let y, — > y 2 — y. The second derivative of the denominator with respect to y\ 
gives 2. Now 

8 2 ,{{8, - 8 2 ) + {y, - y 2 ){8,d 2 - 8 2 2 )) 

= 8, ((dl - d,d 2 ) + (d,d 2 - d 2 2 ) + ( yi - y 2 ){8 2 ,8 2 - 8,81)) 

= d, ((fl? - dl) + ( Vl - y 2 )(8 2 8 2 - 8,81)) 

= (di - 8,81) + (d 2 8 2 - 8,81) + (vi ~ y2)(d!d 2 - 8 2 ,8 2 2 ). 

Evaluating the right-hand side at y — y, — y 2 and noting that 8\8 2 F = 8,8% F at (y,y), 
we just have 8f — 8fd 2 . Hence fl23|) at (y, y) reduces to 

yd% + (c - 1 - y)8 1 8 2 + - d\d 2 ) + ^(28 2 , + 28,8 2 ) - ad, 

= V -(28l + 6dfd 2 ) + (c - 1 - y)8,8 2 + -(25? + 28,8 2 ) - ad,, 
o 4 

where we used 8,F = 8 2 F at (y, y). Comparing the right-hand side with ( 12 ip and by ( 1201) 
we obtain 

|r (y) + (c-i-y) (|r (v) + c -^-f(y) - Y y f{y)) + \ f " {y) ~ P' {y) = °- (24) 

This equation can be written as 

f m (V) = h 2 (y)f"(y) + h,( y )f'(y) + h (y)f(y), 

where 

3(c-l-y) 2 4a 2(c - y)(c - 1 - y) Aa(c-l-y) 

h ^y) = -> h ^y) = — -2 > W = — - 2 — 

are rational functions in y. The coefficient matrix for the Pfaffian system for a one- 
dimensional ODE is simply the companion matrix 

1 
1 

K h (y) h,(y) h 2 (y) / 

Note that the values of /, /', /" and /'" at the origin are given by 

/(0) = F{0, 0) = 1, f (0) = 28,F(0, 0) = -, 

c 

f" (0) = 2a?F(0, 0) + 2 w(o, 0) = |^ + 

3(c) 2 3c(c-|) 

(0) = 2fl?F(0, 0) + 6^ 2 F(0, 0) = 2^ + 6 (^f + ^ ~ t ) ' 

(c) 3 V5(c) 3 5(c) 2 (c-|)/ 

As seen above, the computation using the l'Hopital rule is already tedious for m = 2. 
Actually the computation can be automated by the restriction algorithm for holonomic 
ideals. This will be explained in Section 15.21 
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4 Properties of the Pfaffian system (integrable con- 
nection) for a general dimension 



We now consider our problem for a general dimension. We fully utilize Grobner basis 
theory for the ring of differential operators. In this section we only consider the non- 
diagonal region X . Let K = C(yi, . . . , y m ) be the field of rational functions in y\, . . . , y m 
with complex coefficients. Further let 

R = K(di, ...,d m )= C(j/i, . . . ,y m )(di . . .,d m ) 

be the ring of differential operators with rational function coefficients (see Appendix of 
Nakayama et al. [jjl]). Let / denote the left ideal of R generated by g±, . . . , g m : 

I = (9i, ■ ■ ■ ,9m), (25) 

where is given in 

We now prove the following lemma concerning the commutators of g±, ... , g m . 
Lemma 1. For 1 < i ^ j < m, 

A similar result for 2 F\ is given in Lemma 9.9 of Ibukiyama et al. j?}. Although 
they claim that their Lemma 9.9 follows from a straightforward computation, in fact the 
computation for checking ( 1261) is tedious even for m = 2. However for m = 2, ( 126]) can be 
verified by some software systems (e.g., RisaAsir developing team 0), which can handle 
rings of differential operators. The following program in Risa/Asir 



import ( "names . rr ") ; import ( "yang.rr") ; 
yang.def ine_ring( ["partial" , [yl ,y2] ] ) ; 

Gl=yl*dyl~2+(c-yl)*dyl+(l/2)*(y2/(yl-y2))*(dyl-dy2)-a; 
G2=base_replace(Gl, [[yl,y2] , [y2,yl] , [dyl,dy2] , [dy2,dyl]]) ; 
G=yang.mul(Gl,G2)-yang.mul(G2,Gl)+(l/2)*(yl+y2)/(yl-y2)-2*(Gl-G2) ; 
printf ("G=y,a\n" ,G) ; 



outputs the result G=0. Therefore in the following proof, assuming that (1261) holds for 
m = 2, we show that it holds for m > 2. 

Proof. By symmetry we only need to prove the case % — l,j = 2. Define g±,g2 

gi = y\d\ + (c - y x )d\ + — (#i - d 2 ) - a, 

2 yi - y 2 

~9i = yid\ + (c - y 2 )d 2 + -—^ — (<9 2 - d x ) - a. 

2 2/2 - yi 



Then 



9i=9i + hi, hi = -y^ — — — (di-dk), i = l,2. 
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We already know 

r~ ~1 1 2/1+2/2 ,~ ~s 

l9l ' 92, = -2feT^ (9l - 92) - 

Then 

52] = [01 + ^1, 02 + M = [01, 02] + [/ii, 02] + [01, M + [hi, h 2 ]. 
Therefore it suffices to show 

[^1,02] + [01,^2] + [hi,h 2 ] = -\ , Vl + V \j h 1 - h 2 ). 

2 {yi - y 2 y 

In considering commutators, we only need to look at terms, where a differential oper- 
ator actually differentiate rational functions in yi, . . . , y m . For example consider hig 2 in 
[hi, g 2 ]- In ^102 the only relevant term is di in hi differentiating yi/(yi — y 2 ) in g 2 . Noting 

yi o i 2/2 ,x 2/2 



fli^_ = ^(^!«_-l) 



2/2-2/1 1/2-2/1 (y 2 -y\) 2 ' 

in /i!^ 2 the relevant terms are 

4 (2/2 - 2/1) 2 f^Vi-Vk 4 (2/2 - 2/1) 2 ^ 2/1 - 2/fc 

In (72^1 we need to look at 9i in g 2 differentiating yk/(yi — yk)- Hence we have 

I 2/i Vk 



4 2/2 - 2/i ^ (2/1 - j/a) 
Similarly in [<?i, /12] = — [^2,01] the relevant terms are 

TfL TTl 

4 (2/1 - y 2 ) 2 ^ 2/2 - 2/fc 4 j/ t - ^ (2/2 - 2/fc) 

Finally in [/ii, /12] we look at dk differentiating y^j (yi — yk). Then the relevant terms 

1 Vk 2/2 ,rs o X . 1 Vk 2/1 /o o \ 

4 ^ 2/1 - y* (2/2 - yk) 2 4^?/2- 2/fe (2/1 - 2/*) 

Then the coefficient for —(pi — <9fc)/4 is 

2/2 2/fc ,2/i 2/Jfe 2/i 2/A Vk Vi 



arc 



(2/2 - yi) 2 yi - yk 2/2 - yi (yi - yk) 2 (y\ - y 2 ) 2 y 2 - yk 2/2- Vk(yi- yk) 2 
yi + y2 yk 

(2/2 - 2/i) 2 2/i - 2/fe' 
which coincides with the coefficient of —(d\ — 9fe)/4 in 

1 2/1 + 2/2 , 

2 - 2/ 2 ) 2 

Similarly the coefficients of (d 2 — dk) coincide on both sides. □ 
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We now consider the graded lexicographic term order >-. The initial term of gi (without 
the coefficient yi) is given as 

hv#i = d}. 
We now prove the following theorem. 

Theorem 2. For the term order y, {gi, . . . , g m } is a Grobner basis of I in R and the ini- 
tial ideal is given by (df, . . . , d^) . I is zero- dimensional and the set of standard monomials 
is given by the set of square-free mixed derivatives 

{d h d i2 ...d ik \ 1 < % x < ■ - ■ < i k < m, k < m}, 

which has the cardinality 2 m . 

Proof. By Lemma [T] and the Buchberger's criterion for the ring R (cf. Theorem 1.1.10 
of Saito et al. [23|] ) , gi,i = 1, . . . ,m, form a Grobner basis and the initial ideal is given 
by (df, ...,dl). Let J = (d u ...,d m ). Then J m+1 C (df, ...,8^). Hence / is a zero- 
dimensional ideal. Furthermore this shows that the set of standard monomials is given 
by the set of square-free mixed derivatives. □ 

It follows from Theorem |2] that there exists a Pfaffian system and 2 m x 2 m matrices (as 
Pi(Y) for m = 2 in the expository section [3]) are obtained by the normal form algorithm 
in the ring of differential operators R. The matrices are used to numerically solve the 
associated ODE. However, the derivation of the matrices on computer is heavy and the 
obtained matrices are not in a relevant form for an efficient numerical evaluation. Then, 
we do it by hand in the sequel. 

Consider a higher order derivative . . . d^F of F = i-Fi(a; c;yi, . . . , y m ). If total 
degree of differentiation n = n\ + ■ ■ ■ + n m is greater than or equal to m + 1, then for 
some i we have rij > 2. Then as in the previous section we can use g±F = to decrease 
the total degree of differentiation. Therefore as in (TTTT) . for each m, . . . ,n m , there exist 
2 m rational functions h^ 1 '"^" 1 ^, ij = 0, 1, j = 1, . . . , m, such that 

i i 

flp . . . dfrF = ■ ■ ■ E >''"' <A> ■ ■■ 9 m F - (27) 

ii=0 i m =0 

In the holonomic gradient method, as in the case of m = 2 in ( fl4l) . we only need 
^i™ 1 ' '1'™ where < ni, . . . ,n m < 2 and at most one of n\,...,n m is two, such as 
^h' 1 ' '"' Define a 2 m -dimensional vector of square-free mixed derivatives of F by 
F = (F, diF, d 2 F, did 2 F, . . . , di . . . d m Fy. In F the elements are lexicographically ordered, 
for convenience in programming. diF is written as 

diF = Pi(y)F, i = l,...,m, 

where Pi(y), % = l,...,m, in the Pfaffian system are 2 m x 2 m matrices consisting of 

/4 ni '7 nm) 's. 
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We now study the form of h^ 1 '"^^™ , where < ni, . . . ,n m < 2 and at most one of 
ni, ... , n m is two. Denote [m] = {1, . . . , m}. For a subset J C [m] denote 

dj = U d r 

jeJ 

Choose % G [m] and J C [m] such that i ^ J. Write I = J U {i}. djgiF = djO = 0, where 
gi is in Since i J, we can write djQi as 

ytfdj + (c - yi )dj + ~ V ^(-J^-fl, - QO) - a9j. 



For fc ^ J 



dj{^^(d t -d k )) 



Vk 



Vi - Vk 

On the other hand for k G J, by (J22]) 



JU{fc}, 



'Vi-Vk Vi-Vk 
Here <9j<9fc is not square-free and in fact 



(9/ - <9/<9 fc ) + - 1/8 , 9 (%}uj\W ~ 



(Vi ~ Vk)' 



djd k = did, 



j\{k}, 



which causes recursive application of §!]). In djgi we now separate square- free terms and 
define 

r(i, J;y) — — \{c — y^di - adj + - V yfc (dj - <9/<9 fc ) 



where for J = 0, reflecting the original <7j, we define 



0; y) = - (c - yO^ - a + £ ^ ^ (ft - <9fc 



2 U Vi ~ Vk 



Then dfdjF is expanded as 



= r(i, J; y)F + ^ —3 — {ykdldj\{ k y)F. 



2 t^y*-y k 



(28 



The use of this recursive expression yields an efficient numerical evaluation of the matrices 
of the Pfaffian system. We keep numerical values of d1dj\{k}F in a table and use them to 
evaluate dfdjF and keep it in the table, again. 
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We can also apply the recursion to the last term on the right-hand side. The resulting 
expression for yidfdjF is given as 



dldjF = r(z, J; y)F + lYl ~^—r{K J \ {*}; V)F 

+ i ^ (y-y k L _ yk , r(h,J\{k 1 ,k 2 };y)F 



fci , fco :distinct 



1 1 

a E 7 w 77 T r ( fc 3, -/\ {fci,^2,fc 3 };2/)F+ . 



fc^ , &2 ,^3 : distinct 



+ 2 ' J| fcl ,5, 6J (vi-yki)(yki-Vk 2 ) ■■■(yk ]J] ^-yk m ) r ^ klJ],(l),v ^ R ^ 

fcj^ , .. .,fci j- 1 : distinct 

Now in (@J we write Yr l jl = (3 = . . . ,/3 m ), where . . . , (3 m are distinct, and 
define a 2 m -dimensional vector valued function G in a scalar x by 

m 

= exp(-x^ f3i)x mn/2 F{Px). 

i=l 

Then G satisfies the ODE 

^ = HE a) 7 *- + ^ J 2- + E wa) £ ( 30 ) 

where is the 2 m x 2 m identity matrix. We denote the right-hand side as PpG. We 
now prove the following theorem, which is important for guaranteeing stability of ODE 
at x = +oo. 

Theorem 3. As x — > oo 

Pp = A + O(l/x), 

where A only depends on (3 and the 2 m eigenvalues of A are given as —t\{3\ — - e m (3 m , 

where (d, . . . , e m ) E {0, l} m . 

Proof. Note that y\ — fliX, . . . , y m — f3 m x = 0(x). Divide fl29l) by yi = fyx. Then on the 
right-hand side of (]29|) . the only constant order term is d[ in r(i, J; y). Now 



A m 

—d I F((3x) = J2^d I F(Px) 



dx 

i=l 



E Pid*d im F(Px) + ]T Pid IU{ i } F(l3x) 
(J2^)diF((3x) + 0(l/x) +J2Pi d iu{i}F(Px). 



iei i&I 
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This implies that the 7-th diagonal element of Aq is given 

m 

-I>+Ea = -£a- 

i=l iel i0 

Furthermore the (I, IU {z})-element of Aq is Other elements of Aq are zeros. Hence Aq 
is an upper triangular matrix with diagonal elements I ^ [m]. The theorem 

holds because the diagonal elements of an upper triangular matrix are its eigenvalues. □ 



5 Some results of symbolic computation 

In this section we present some results on symbolic computation for the initial values (cf. 
(|17p ) and the restriction for diagonal regions (cf. Section 1X2]) . We omit writing down the 
fully expanded form of (|29|) . since the recursive formula (1281) can be directly used in our 
implementation of holonomic gradient method. 



5.1 Initial values 

Initial values for our holonomic gradient method can be obtained by expressing iF\ in 
terms of monomial symmetric polynomials as in ([16]) . We denote the relation between 
the zonal polynomials and the monomial symmetric polynomials in ( lT5"j) as 



C K (Y) = J2^,xM x (Y) 1 



where A <j k means that A is dominated by k, i.e. Yli=i — ^2i=i K i f° r a ^ s - Then iFi 
is expressed as 



(c) K k\ ' 

k=0 Xhk Khk,K>\ v > K 



1 F 1 (a;c;F) = ^^? A (a,c)M A (y), qx(a,c) = ^ 



A recurrence relation for c Kj a's is given by James [9| (see also (14) in Section 7.2.1 
of Muirhead 18j] and Section 4.5.4 of Takemura 30|), which can be used to compute 
q K (a, c). However James' recurrence relation works for each C K separately. Recently Koev 
and Edelman 10| gave a much improved algorithm based on recursive relations among the 
values of zonal polynomials for m variables and m — 1 variables. For our implementation 
of holonomic gradient method, we adapted Koev-Edelman's recurrence relation also for 
derivatives of ±Fi to evaluate the initial values. 

Close to the origin, we can use rough initial values given by the linear approximation 
as in (fTSI) . Then we only need k = (ki, ■ ■ ■ ,h) such that ki = ■■■ = ki = 1 or kx = 2, 
k\ = ■ ■ ■ = k\ = 1. Some q\(a, c)'s for small A's are as follows. 

, a (a) 2 (a) 2 , 2(a) (ljl) (a) 3 2(a) (2jl) 

(a) 3 3(o)(2,i) (a)(i,i,i) (Q)4 4(a) (2 ,2) H(a)(3,i) 2(a) (2 ,i,i) 

15(c) 8 + 5(c) (2)1) + 3( C ) (1>ltl) ' g(2 ' M) 70(c) 4 + 45(c) (2)2) + 63(c) (3)1) + 9(c) (2>1)1) 
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where stands for the unique partition of zero. Write (l k ) = (1, . . . , 1), (2, l fc ~ 2 ) = 
(2, 1, . . . , 1), which are partitions of k. Given the above quantities, the linear approxi- 
mation of di . . . d[F(Y), < I < m, for Y — (y±, . . . ,y m ) close to the origin is expressed 
as 

<9i . ..diFiY) = (a, c)+2g (2> ii-i) (a, c)(j/H \-yi) + ( a , c ) (yi+i H hy m ), (31) 

where for / = the second term on the right-hand side is zero and for I = m the third 
term is zero. We found that initial values b y (13 ip are practical enough for m < 5. 

In fact, by Lemma 1 in Section 4.5.2 of 3o[ and by Proposition 7.3 of 25|, qnk\(a, c) 
and 5(2 ifc-2)(a, c) are explicitly written as follows: 



q (1 k)(a,c) = Tk\ > = J ~\ )-$-, 

g (2 lfc - 2) (a, c) = 2 fc (fc - 2)! V O^g^ I 2 ^ I ! + j ( f f] + ^ ^ _ <)) 



where is the length (number of non-zero parts) of k — (hi, ... , 

For larger values of m we need higher order terms for initial values. For two partitions 
/i, A, we write fi C A to denote /ij < Aj for all z. For two partitions k, is, we denote by 
Kl±ly the concatenation of k and z/ obtained from (k±, v\, k 2 , z/ 2 , . . .) by sorting. Consider 
a rectangular partition r = (t, . . . , t) = (t l ) h £/. For r = (£*) and A D r we define 

/(A, r) = {(k, is) I k ttl v — A, r C k, = 0}. 

Consider d"M\(Y) = ■ ■■d£>Mx(Y). Note that d"M\(Y) = 0, if // £ A. For a 

rectangular r = (£') we can calculate 9 t A^a(^) by the following lemma. 

Lemma 2. For r = (t l ) C A 



d T M x (Y)= V — M K - T (y yi)M u (yi +1 ,...,y m ), 

(k,i/)G/(A,t) V ' 



where 7! = H^TiO / or a partition 7, and n — t = («i — t, . . . , ki — t) . 

Proof is straightforward and omitted. Using this lemma we have the following propo- 
sition. 



Proposition 1. For a rectangular partition r = (t l ), 

7! 

( 7 -r)! 



9 r iF!(a;c;F) = J^^g A (a,c) ^' M^ T (y u . . . , yi)M v (yi+\, ■ ■ ■ , y m ). 



k=tl Ahfc, (7,!/)eT(A,r) 
rCA 



(32) 
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For our initial values we only need to consider r = (1 ). We obtain fl3Tl) if we only look 
at linear terms on the right-hand side of fl32l) . Note that since i-Fi(a; c; Y) is a symmetric 
function in yi, . . . , y m , other derivatives are obtained by permutation of yi, . . . , y m . 

Although f[3"2"l) only gives derivative with respect to a rectangular partition r = (t l ), we 
can obtain other derivatives d 1 ^ 1 . . . dj J ' h iF 1 (a; c;Y), jix > ■ ■ ■ > /i h , by repeated application 
of (1321) for different values of Vs. 



5.2 Restriction to diagonal regions 



As mentioned at the end of Section [3721 the tedious operation involving the l'Hopital rule 
for the diagonal region can be performed by the restriction algorithm for holonomic ideals. 
The following program in Risa/Asir for m = 2 



import ("names .rr") ; import ("nk_restrict ion . rr" ) ; 

Gl=yl*dyl~2 + (c-yl) *dyl+(l/2) * (y2/ (yl-y2) ) * (dyl-dy2) -a; Gl=red( (yl-y2) *G1) ; 
G2=base_replace(Gl, [[yl,y2] , [y2,yl] , [dyl,dy2] , [dy2,dyl] ] ) ; 
F=base_replace ( [Gl , G2] , [ [y 1 , y] , [y2 , y+z2] , [dy 1 , dy-dz2] , [dy2 , dz2] ] ) ; 
A=nk_restriction.restriction_ideal(F, [z2,y] , [dz2,dy] , [1,0] I param= [a, c] ) ; 



produces the output 

-y"2*dy"3+(3*y"2+(-3*c+l)*y)*dy~2+(-2*y"2+(4*a+4*c-2)*y-2*c"2+2*c)*dy-4*a*y+(4*c-4)*a 

This is the same as (124|) . Adapting the above program for m = 3, we obtain 

y 3 f""(y) + (-6y 3 + (6c-4)y 2 )r(y) 

+ (Hi/ 3 + (-10a - 22c + lS)y 2 + (11c 2 - 17c + 4)y)f"{y) 

+ (-6y 3 + (30a + 18c - 18)y 2 + ((-30c + 34)a - 18c 2 + 34c - 12)y 

+ 6c 3 - 16c 2 + 10c)f'{y) 
+ (-18ay 2 + (9a 2 + (36c - hl)a)y + (-18c 2 + 48c - 30)a)/(y) = 0. 

For m = 4, we found that the computation by Risa/Asir takes too much time and 
memory. 

We conjecture that the ideal / generated by n^i^ ~ Vi)9h i = in the 

Weyl algebra D m is an holonomic ideal. In fact, the conjecture can be checked for small 
dimensions m on a computer. See the Appendix fAl If I is a holonomic ideal, then 

J=(i + (vi- 2/2) An + (yi - 2/3) An h 1- (yi - y m )D m ) n C(y 1} d yi ) 

is not and is an holonomic ideal in A by the theorem of Bernstein (see, e.g., the 



Chapter 5 of 23]). The generators of J is ordinary differential equations for the function 
restricted to the diagonal y\ — ■ ■ ■ — y m . Thus, the holonomicity implies the existence 
of the diagonal ordinary differential equation. The ideal J can be obtained by Oaku's 
algorithm (jioj) based on Grobner bases and the Risa/Asir package nk_restriction 
uses this algorithm. 
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6 Numerical experiments 



We implemented the holonomic gradient method in a straightforward manner. Our source 
programs in the language C are available from 



http: //www. math. kobe-u.ac . jp/OpenXM/Math/lFl 



The updating step of the holonomic gradient method was implemented using the 
recursive relation ( 128]) for a general dimension. For initial values we adapted Koev and 
Edelman 1(1 for derivatives of \Fi as discussed in Section 15.11 



The accuracy of the holonomic gradient method can be simply checked by looking at 
the numerical convergence Pr[£i < x] — > 1 as x — > oo. This is because the initial values 
are evaluated at small x > and Pt[£i < x] at large x is obtained after many updating 
steps. This is another advantage of our method. 

Also we can use the following simple bounds for the upper tail probability for the 
purpose of checking. Let Pt[£i < denote the probability under the covariance matrix 
S. Consider £ = diag(<r 2 , . . . , cr 2 J, of > • • • > cr^. Then by standard stochastic ordering 
consideration, we have 

Pr[£t < x\ diag(a 2 , . . . , af)} < Pr[£i < x\ diag((x 2 , . . . , <j 2 m )\ 

< Pr[li < x\ diag((T 2 , 0, . . . , 0)]. (33) 

The upper bound coincides with the cumulative probability of chi-square distribution 
with n degrees of freedom (cf., j28|,|3lj|). Accurate approximation for the lower bound 
Pr[£i < x\afl m } is given by the tube method ([HI], 0)- 

We first consider the case m = 5,n = 7, S~ x /2 = (3 = (1,2,3,4,5). For x = 20, 
the two bounds in ( 13"3"|) are given as 0.9996034 and 0.9999987. With the initial value of 
Xq = 0.01 and step size 0.0001 we obtained 

Pr[£i < 20] = 0.999972. 

The cumulative distribution function for this case is plotted on the left part of Figure [TJ 
Next we consider the case m — 10, n — 12, /3 = (1, 2, . . . , 10). For x = 30, the bounds 
are (0.99866943,0.99999998). For generation of initial values it takes about 20 seconds 
for approximating iFi and its derivatives up to the degree 20 with an Intel Core i7 CPU. 
With the initial value of xo = 0.2 and step size 0.001, we obtain 

Pr[£x < 30] = 0.999545 

in about 75 seconds. The cumulative distribution function for this case is plotted on the 
right part of Figure [TJ We see that enough accuracy is obtained even for m = 10 within 
practical amount of time. 

The complexity of numerically solving the ODE for G ( 130]) is 

0(m2 m ) x (steps of the Runge-Kutta method with a prescribed precision). 

In fact, since the matrix Pi(j3x) has sparsity, each vector Pi(j3x)G(x), which has 2 m 
elements, can be evaluated in 0{2 m ) steps at x from the values of G(x) by utilizing 
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Figure 1: Cumulative distributions for m = 5 and m = 10 



7 Discussion of open problems 



The holonomic gradient method ([19J) gives a general algorithm for obtaining the partial 
differential equations satisfied by parametrized definite integrals such as the normalizing 
constant of a family of probability distributions. In fact, in 



19 and Sei et al. 24 



we 

used the holonomic gradient method for deriving the partial differential equations of 
the normalizing constants and for maximum likelihood estimation for distributions in 
directional statistics. For the case of \Fi, the partial differential equations were already 
derived by Muirhead 16J more than 40 years ago. Our use of those partial differential 
equations for numerical evaluation of is very straightforward as discussed in Section [3] 
for the two dimensional case. Yet, from the viewpoint of holonomic functions, the partial 
differential equations of Muirhead 16j present many interesting open problems. 

One important question is to obtain the ordinary and partial differential equations 
for the diagonal case as discussed in Section 13.21 for the case of m — 2. For a general 
dimension m > 2, it is desirable to be able to handle various patterns of diagonalization, 
such as the two-block diagonalization yi = ■ ■ ■ = yi > yi+i = ■ ■ ■ = y m . A direct "by 
hand" calculation using the l'Hopital rule becomes quickly infeasible when we increase m. 
Also the use of the restriction algorithm for holonomic ideals is limited by computational 
complexity. It is in fact a very heavy algorithm. Currently the nk_restriction routine of 
Risa/Asir in Section I5T21 takes too much time for m > 4. One possibility is to follow the 
approach in Muirhead 16| and Sugiyama et al. 29j, where differentiation with respect to 
elementary symmetric functions of the roots of Y are considered. As discussed in Section 
15. 2\ we conjecture that the ideal I generated by Ylj^ilJi ~~ Vj)9h i — l,---,m in the 
Weyl algebra D m is an holonomic ideal. Holonomicity guarantees the existence of partial 
differential equations for diagonal regions. 

Another question is to consider the asymptotics for Pr[£i > x] = 1 — Pt[£i < x] as 
x — > oo. As mentioned in the previous section, this tail probability can be approximated 
by the tube method ([ll|, 12j])- O ne theoretical problem in applying the tube method 
is that only the approximation for the tail probability itself has been justified and the 
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justification of its derivatives has to be proved. However it is obvious that the current 
approach of taking the initial values close to the origin causes difficulty in precision for 
the extreme upper tail probability, in the case we want to evaluate the small probability 
Pr[£i > x]. Hence it is desirable to be able to use tube formula approximation as the 
initial values at x = oo. 

From computational viewpoint, our holonomic gradient method has exponential com- 
plexity in the dimension m. We need to keep the 2 m -dimensional numerical vector F 
in memory at each step of the iteration. For m = 20, the dimension of the vector is 
about one million. Hence we do not expect that the current implementation of the holo- 
nomic gradient method works for m = 20. It might be possible to improve our current 
implementation by fully exploiting the fact that \F\ is a symmetric function in Y. 



A Holonomicity for dimension two 

In the theory of holonomic functions, the holonomicity of the left ideal generated by the 
set of partial differential operators is an important question. In fact, the existence of the 
ordinary differential equation with polynomial coefficients for the function restricted to 
the diagonal region follows from the holonomicity. Holonomicity of the ideal generated by 
c/i, g 2 in the two-dimensional case can be verified by Grobner basis computation. Here we 
present this result. As to a general introduction to holonomic ideals and Grobner bases, 



we refer to the Chapter 1 of Saito et al. [23 



Note that the holonomicity on the non-diagonal region X follows from Theorem [2] 
because the zero set of yiQ — 0, % — 1, . . . ,m contains the characteristic variety on X. 
The holonomicity on X can also be proved by an analogous method with Ibukiyama et al. 

0. 

Let D 2 be the second Weyl algebra. For P = J2 k =o E Ql+Q2 =fc fa uaa (yi, V2)d^d 2 2 G 
D 2 , we define in(P) by 

in(P) = /«!,« 2 (!/i.^rS 2 e%.2/2,(i,(2], 

ai+ot2=d 

where we assume that f ai! a 2 (yi, y 2 ) € C [2/1,2/2] an d that / Ql ,« 2 (2/i, 2/2) 7^ for some cti,a 2 
with ol\ + a 2 = d. For a left ideal / of D 2 , the characteristic variety ch(7) is defined by 

ch(/) = {(2/1,2/2,6,6) e C 2 - 2 1 VP e I, in(P)( yi , 2/ 2 , 6, 6) = 0}. 

It is a basic fact that the dimension of the characteristic variety ch(J) of the proper left 
ideal / of P 2 is greater than or equal to 2. A left ideal / of D 2 is called holonomic if the 
dimension of the characteristic variety ch(J) equals 2. 

Let Pi = (2/1 — 2/2)2/1, P2 = (2/2 — 2/1)^2 and let / be the ideal of D 2 generated by Pi and 
P 2 . We will show that / is holonomic. Let S = y 2 d 2 P 1 + y\d\P 2 + c[d 2 P x + d 1 P 2 ) G I. By 
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direct calculation we have 

S =(yh2 - Vivl + f - 2y m )dld 2 + {-yly 2 + yi y\ - 2 ym + V l)d x d\ - V \d\ - |«9 2 3 
+ (ayf - ay x y 2 - - yi)d\ + {-y 2 + 2ay 2 - ay x y 2 + ay\ - -y 2 -)^ 

2 2 3c?/! 3«/ 2 

+ + 2cyi?/ 2 - cy 2 + Ay x y 2 - — + y x + y 2 )di d 2 

+ (aq/i — aq/ 2 + 2cw/i + q/i — c 2 )d\ + (— c 2 — aq/i + acy 2 + cy 2 )d 2 + 2ac. 

Hence 

in(5) =(y 2 2/ 2 - ywl + f - 2j/ lWa )£6 + H&fa + Vivl ~ 2 Vl y 2 + |%£ 2 2 - y£i 3 - f £ 2 3 - 

Consider the ideal J of C[y±, y 2 , £1, £ 2 ] generated by in(Pi), in(P 2 ) and in(5). The following 
is a Grobner base of J with respect to the graded reverse lexicographic order with yi > 
Vl > Ci > 6: 

{2/ 2 £i£ 2 + Hild + 3y 2 2 £i£ 2 + y 2 2 e 2 , ymg + 3?/iy 2 £ 2 6 + + y 2 2 & 

Thus the Krull dimension of J is 2. Since {Pi, P 2 , S} C /, the characteristic variety ch(J) 
of / is contained in the zero set of J. This implies that the dimension of ch(J) does not 
exceed 2 and hence the dimension of ch(J) is equal to 2. Therefore / is holonomic for 
m = 2. 



B R source program for dimension two 

The following program in data analysis system R implements holonomic gradient method 
for m = 2 of Section 13.11 based on deSolve add-on package for R. 

library (deSolve ) 

m <- 2 # dimension 

n <- 3 # degrees of freedom 

x <- 4.31600 # specify x. We evaluate Pr( 11 < x ) 
bl <- 1; b2 <- 2 # (bl,b2) = (1/2) diag(Sigma"{-l}) 

a <- (m+l)/2; c <- (n+m+l)/2; totalsteps <- 10000; stepsize <- x/totalsteps 
# h's 

h2000 <- function(yl,y2) a/yl 

h2010 <- function(yl,y2) -(c-yl)/yl - y2/ (2*yl* (yl-y2) ) 
h2001 <- function(yl,y2) y2/ (2*yl* (yl-y2) ) 
hl200 <- function(yl,y2) a/ (2*y2* (y2-yl) ) 

hl210 <- functional, y2) 3/ (4* (y2-yl) '2) + a/y2 - (c-yl) / (2*y2* (y2-yl) ) 

hl201 <- functional, y2) -3/ (4* (y2-yl) ~2) 

hl211 <- functional, y2) -(c-y2)/y2 - yl/ (2*y2* (y2-yl) ) 

#initial values 

xl <- bl*stepsize; x2 <- b2*stepsize 

fi <- c(a/c, (a*(a+l))/(c*(c+D) , (a* (a+1) ) / (3*c* (c+1) ) + (2*a* (a-1/2) ) /(3*c* (c-1/2) ) , 
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(a*(a+l)*(a+2))/(5*c*(c+l)*(c+2)) + (4*a* (a+1) *(a-l/2) ) / (5*c* (c+1) * (c-2/1) ) ) 
fi <- c(l+(xl+x2)*fi[l] , fi[l]+xl*fi[2]+x2*fi[3] , f i [1] +x2*fi [2] +xl*f i [3] , f i [3] + (xl+x2) *f i [4] ) 



# gradient 

fllm2 <- function(y,fv,parm){yl <- y*bl ; y2 <- y*b2; 
list(c( 

bl*fv[2] + b2*fv[3] , 

bl* (f v [1] *h2000 (yl ,y2) +f v [2] *h2010 (y 1 , y2) +f v [3] *h2001 (y 1 , y2) ) + b2*f v [4] , 
b2* (f v [1] *h2000 (y2 ,yl) +f v [2] *h2001 (y2 ,yl) +f v [3] *h2010 (y2 , y 1) ) + bl*f v [4] , 

bl*(fv[l]*hl200(y2,yl)+fv[2] *hl201(y2,yl)+fv[3] *hl210(y2,yl)+fv[4] *hl211 (y2 ,yl) ) 
+b2*(fv[l]*hl200(yl,y2)+fv[2] *hl210(yl ,y2)+f v [3] *hl201 (yl ,y2) +f v [4] *hl211 (yl ,y2) ) ) 

)} 

output <- ode (f 1 ,f unc=f llm2 , (1 : totalsteps) *x/totalsteps) 

probO <- ((bl*b2)"(n/2)*gamma(a)*gamma(a-l/2))/(gamma(c)*gamma(c-l/2)) * x~(n*m/2) * exp(-x* (bl+b2) ) 
cat("x=",x, "prob=", output [totalsteps, 2] *probO, "\n") 
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